# 消费和储蓄决策
rm(list = ls())
library(pacman)
p_load(bvpSolve,ggplot2,magrittr,signal,rootSolve,matlab)
devtools::load_all()

# 参数设置
beta <- .97
R <-  .02
WL <- 10
numT <- 31 # 31期

find_c <- function(cc, parms = list(numT = 30, beta = 0.97, R = 0.02, WL = 10)){
  eps <- B <- numeric(parms$numT)
  # 初值：B0 = 0
  B[1] <- parms$WL - cc[1]
  # B迭代
  for (i in 2:parms$numT) {
    B[i] <- (1+parms$R)*B[i-1] + parms$WL - cc[i]
  }

  # 消费迭代
  for (i in 1:(parms$numT - 1)) {
    eps[i] <- cc[i+1] - parms$beta * (1+parms$R) * cc[i]
  }
  # 终值：B(T) = 0
  eps[parms$numT] <- (1+R)*B[parms$numT-1] + parms$WL - cc[parms$numT]
  return(eps)
}

# 解多元非线性方程组
cc <- 10*rep(1,numT)
ans <- multiroot(find_c, start = cc, parms = list(numT = numT, beta = beta, R = R, WL = WL))
picdata <- data.frame(ts = 1:length(ans$root), cc = ans$root)
ggplot(picdata, aes(x = ts, y = cc)) + geom_line()


# 稳态
A <- matrix(c(beta*(1+R),0,-beta*(1+R),1+R),nrow = 2, byrow = T) - eye(2)
eigen(A)
-solve(A) %*% matrix(c(0,WL),ncol = 1)

# 这本质上是一个边值问题，因此可以直接用边值问题来解
cns_sv <- function(x, y, pars){
  # y[1]:c
  # y[2]:B
  list(c(beta*(1+R)*y[1]-y[1],
         (1+R)*y[2]+WL-y[1]-y[2]))
}
ans <- bvptwp(yini = c(NA,0),yend = c(NA,0),func = cns_sv, x = seq(1,32))
